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We present a continuous formulation of epidemic spreading on multilayer networks using a tenso- 
rial representation, extending the models of monoplex networks to this context. We derive analytical 
expressions for the epidemic threshold of the SIS and SIR dynamics, as well as upper and lower 
bounds for the disease prevalence In the steady state for the SIS scenario. Using the quasi-stationary 
state method we numerically show the existence of disease localization and the emergence of two or 
more susceptibility peaks, which are characterized analytically and numerically through the inverse 
participation ratio. Furthermore, when mapping the critical dynamics to an eigenvalue problem, we 
observe a characteristic transition in the eigenvalue spectra of the supra-contact tensor as a function 
of the ratio of two spreading rates: if the rate at which the disease spreads within a layer is compara¬ 
ble to the spreading rate across layers, the individual spectra of each layer merge with the coupling 
between layers. Finally, we verihed the barrier effect, i.e., for three-layer conhguration, when the 
layer with the largest eigenvalue is located at the center of the line, it can effectively act as a barrier 
to the disease. The formalism introduced here provides a unifying mathematical approach to disease 
contagion in multiplex systems opening new possibilities for the study of spreading processes. 


I. INTRODUCTION 

Epidemic like spreading processes are paradigmatic, as 
they can describe not only the temporal unfolding and 
evolution of diseases, but also of ideas, information and 
rumors in fields as diverse as biological, information and 
social sciences [T]. Due to their fundamental nature and 
simplicity, two particular models have received special 
attention by the scientific community, the susceptible- 
infected-susceptible (SIS) and the susceptible-infected- 
recovered (SIR). In both models, an infected individual 
spreads the disease to its neighbors at a given (spread¬ 
ing) rate and infected individuals recover at some other 
rate. The difference between both scenarios lies in the 
fact that in the SIS case, once recovered, infected individ¬ 
uals can catch the disease again, and, therefore, they go 
back to the susceptible state. On the contrary, in the SIR 
model, recovered individuals are supposed to acquire per¬ 
manent immunity and do not play any active role in the 
spreading process anymore. There are many other vari¬ 
ations of these two models, including more realistic and 
intricate compartmental models [T]. However, these two 
schemes are sufficient to capture the main phenomenol¬ 
ogy of disease dynamics — and many other contagion 
like processes — including the onset of epidemics, while 
remaining simple. 
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Originally, the modeling of diseases was confined to ho¬ 
mogeneous systems, where any pair of individuals have 
the same contact probability EE]. However, most real- 
world networks are heterogeneously organized, which led 
to reexamine previous results considering non-trivial pat¬ 
terns among individuals, such as power-law degree distri¬ 
butions E^6]. In Cl, the authors presented the hetero¬ 
geneous mean-field approach (HMF), showing that the 
epidemic threshold tends to zero in the thermodynamic 
limit on scale-free networks when they characteristic ex¬ 
ponent is less than 3. This observation about the role 
of network organization changed completely our previ¬ 
ous understanding of how disease outbreaks should be 
modeled and controlled, placing the focus of attention 
not only into new ways to model disease dynamics, but 
also into the incorporation of real contact patterns in the 
dynamical settings EIHHII]. 

Since then, many computational and theoretical frame¬ 
works have been proposed, which undoubtedly had made 
the modeling of disease contagion an active area of re¬ 
search and have provided new phenomenological insights 
and accurate methods for the study of real outbreaks. 
For instance, instead of the HMF approach, one can 
adopt the quenched mean field (QMF) method, where 
a specific network is fixed and the dynamics is modeled 
in terms of nodal probabilities min]. The results ob¬ 
tained with the latter approach show that the epidemic 
threshold depends on the inverse of the leading eigen¬ 
value of the adjacency matrix m [13] — a similar re¬ 
sult was also obtained using a discrete Markov chain ap- 
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proach El- Other scenarios explored recently include 
the case of temporal networks US HU, competing and 
interacting diseases [HHIS] as well as the inclusion of hu¬ 
man behavioral responses pH - ES] , 

However, the vast majority of the works so far deal 
with single-layered networks, despite the fact that many 
real systems exhibit a large degree of interconnectivity 
and hence should be modeled as multilayer networks m 
Such systems represent multimodal, multicategorical or 
temporal interactions, as for instance social relations, the 
ecosystem formed by different online social networks or 
modern transportation systems m- Cozzo et al. [25] 
showed that disregarding the multilayer structure can 
lead to misleading conclusions, missing fundamental as¬ 
pects of the critical dynamics of spreading-like processes. 
Such findings reinforce the importance of a more detailed 
investigation of contagion processes on multilayer net¬ 
works. Here, we develop a theoretical and computational 
framework for the analysis of disease spreading, general¬ 
izing the results of Ref. m to multilayer networks. A 
continuous counterpart to the model presented in |28| is 
provided in terms of the tensorial notation introduced 
in |23]. Our methodology allows for several new results. 
First, we are able to write down in a compact form the 
equations describing the disease dynamics in a multilayer 
system. Secondly, we derive the corresponding epidemic 
thresholds for the SIS and SIR cases as well as estab¬ 
lish bounds for the prevalence of the disease in the SIS 
scenario. We also identify previously unnoticed multiple 
susceptibility peaks and disease localization, which are 
traced back to the very topological nature of the system 
and described in terms of the eigenvalue spectra of the 
supra-contact tensor and the localization of eigenstates. 

The rest of the paper is organized as follows: we first 
formally define the concept of multilayer network, intro¬ 
ducing the tensorial notation. Next we derive the equa¬ 
tions describing the dynamics of the disease for the SIS 
scheme, calculating the upper and lower bounds for the 
prevalence of the disease in the steady state, followed 
by the analytical expression for the epidemic threshold, 
which is also derived for the SIR model. Furthermore, 
we use the results in m to define some constraints 
on the critical point. In addition, we explore the no¬ 
tion of localization of eigenstates, formerly applied on 
epidemic spreading in m, to inspect localization tran¬ 
sitions, which were verified by multiple susceptibility 
peaks. Finally, we also present results from extensive nu¬ 
merical simulations considering multiplex networks with 
scale-free and scale-rich structures, computing their re¬ 
spective epidemic thresholds. Finally, we present our 
conclusions in the last section. 


II. CONTINUOUS FORMULATION FOR 
MULTILAYER EPIDEMIC SPREADING 

Multilayer networks have been shown to better de¬ 
scribe interdependent systems. Mathematically, they can 


be described by either generalizing the matrix represen¬ 
tation and formalism or by encoding the system’s 
topology in a tensorial representation, which was recently 
proposed [23] and first applied to describe a dynamical 
process in [32]. Here, we use the latter framework to 
formulate a continuous time Markov chain model that 
describes the evolution of an epidemic processes. 

A. Tensorial representation 

Tensors are elegant mathematical objects that gener¬ 
alize the concepts of scalars, vectors and matrices. A 
tensorial representation provides a natural and concise 
framework for modeling and solving multidimensional 
problems and is widely used in different fields, from lin¬ 
ear algebra to physics. In particular, general relativity is 
completely formulated under the tensorial notation. Here 
we use the representation formerly presented in [22]. We 
also adopt the Einstein summation convention, in order 
to have more compact equations: if two indices are re¬ 
peated, where one is a superscript and the other a sub¬ 
script, then such operation implies a summation. Aside 
from that, the result is a tensor whose rank lowers by 
2. For instance, A^A'^ = notation we 

use greek letters to indicate the components of a tensor. 
In addition, we use tilde (~) to denote the components 
related to the layers, with dimension m, while the com¬ 
ponents without tilde have dimension n and are related 
to the nodes. 

A multilayer network is represented as the fourth-order 
adjacency tensor M G ]^"xnxmxm^ which can represent 
several relations between nodes [22] 

m 

= E C<^(hk)E\(hk) = 

m n ^ ' 

= E E Wij{hk)£^^{ijhk), 

h,k=i 

where Ej{hk) G g j^nxnxmxm 

indicate the tensor in its respective canonical basis. Ob¬ 
serve that we can extract one layer by projecting the 
tensor to the canonical tensor Ej{ff). Formally, 

from |22] we have 

M^^Ejirf) = C^“(rT) = A^(r), (2) 

where r G {1,2, is the selected layer and A^{f) is 

the adjacency matrix (rank-2 tensor). Moreover, aiming 
at having more compact and clear equations we define 
the all-one tensors G M" and G Here, we 

restrict our analysis to multilayer networks with a diago¬ 
nal coupling In other words, each node can have at 
most one counterpart on the other layers. In addition, for 
simplicity, we focus on unweighted and undirected con¬ 
nected networks, in which there is a path from each node 
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to all other nodes. For complementary information about 
the tensorial representation, its projections and the gen¬ 
eralization of the eigenvalue problem, see Appendix [X| 


B. The Susceptible-Infected-Susceptible (SIS) 
model 


Despite its simplicity, the susceptible-infected- 
susceptible (SIS) and susceptible-infected-recovered 
(SIR) models capture the main features of disease 
spreading [T]. In this section we focus on the first 
order approximation of the SIS model. Additionally, we 
present some aspects of the SIS exact formulation on 
Appendix |B 1| and a brief analysis of the SIR model on 
Appendix |C[ 

We model the SIS disease dynamics associating a Pois¬ 
son process to each of the elementary dynamical transi¬ 
tions: intra and inter layer spreading and the recovery 
from the infected state. The first two processes are as¬ 
sociated to the edges of the graph and are characterized 
by the parameters A and rj, respectively. The latter tran¬ 
sition is modeled in the node, also via a Poisson process 
with parameter S. Using the tensorial notation defined 
above, the equations describing the systems dynamics 
read as 


dX 


ps 


dt 




ps 


(> - 


A7^“|(A,^?)A„y, 


(3) 


where the supra contact tensor is defined as 

- S]), 

(4) 

which encodes the contacts. It has a similar role as the 
matrix R in [5S] . Notice that we have implicitly assumed 
that the random variables X^g are independent. For¬ 
mally, if the state variable (Bernoulli random variable) 
Sps is such that Sp-s = 1 when the node /3 on layer 5 is a 
spreader and S^g = 0 otherwise, then P[Spg = 1] = 

In this way, the independence of random variables implies 
that P[Sfgg = = 1] = = 1]P[S^^ = 1] = 

X^gX^p ■ Gator and Van Mieghem [33] proved rigorously 
that the states of any two nodes in the SIS model are 
non-negatively correlated for all finite graphs. This re¬ 
sult can be easily extended to our case, since we are con¬ 
sidering constant rates and Markovian processes. Due to 
the positive contribution of the infected nodes we have 
P[Spg = = 1] > P[Sfgg = 1], implying that the 

model is always overestimated. A similar conclusion was 
also obtained in m for the monolayer case. 

Naturally, the order parameter, also called macro-state 
variable, is defined as the average of the individual prob¬ 
abilities, formally given by 


p = —XpgU^~^. 

nm ^ 


(5) 


Note that the steady state is not an absorbing state in the 
Markov sense, since there is a set of possible states where 


the system remains trapped and there is a stochastic vari¬ 
ation over time. In addition, note that there are many 
different configurations for which the fraction of infected 
nodes is the same. More formally, there is a set of states 
above the threshold, which have finite probability larger 
than zero, configuring a meta-state. The only absorbing 
state of this set of equations is thus the disease-free state, 
since when it is reached the (micro and macro) dynamics 
stops. 

Furthermore, one of the most important concepts on 
disease spreading processes is the epidemic threshold: be¬ 
fore the threshold, the system is in a disease-free state. 
On the other hand, when increasing the spreading rate 
it drives the population to an endemic state. In other 
words, there is a nonzero probability that the disease re¬ 
mains on the population, configuring the meta-state de¬ 
scribed above. Analogously to the results for monolayer 
systems we have a critical point given as 



where Ai is the largest eigenvalue of TZ. The com¬ 
plete derivation of the critical point is presented in Ap¬ 
pendix |B 2| Observe that the eigen-structure of the ten¬ 
sor TZ is the same as for the matrix R in |3S], since 
it can be understood as a flattened version of the ten¬ 
sor 7?,“J(A,?7). As argued in jH], the supra-adjacency 
matrix corresponds to a unique unfolding of the fourth- 
order tensor TZ yielding square matrices. Moreover, if 

r]M^TE^{j3(3) <C XM^^Ej{SS), the critical point is dom¬ 
inated by the individual layer behavior and the epidemic 
threshold is approximated to that of a SIS model on 
monolayers, when considering the union of m disjoint 
networks. Consequently, the epidemic threshold is deter¬ 
mined by the largest eigenvalue, considering all layers. 
The same conclusion was reached in |28| using perturba¬ 
tion theory on the supra-contact matrix. 

Finally, the nodal probability on the steady state can 
be bounded by 


1 - 


1 




[(J) 


A ^min _ I 


<X^g<l- 


-]dps + l 


ii) 


(7) 

where denotes the probability that node /3 in layer 

6 is in the steady state regime, = TZ°‘^(X,r])Uap (also 

defined in B9) and d™'" = Min{d^^}. The derivation of 
such bounds are shown in details on Appendix |B 3| In¬ 
terestingly, observe that the higher d“™, the closer the 

lower and upper bounds. In the extreme case ^^^—>-00 


the bounds approach each other and all nodes tend to be 
infected. Phenomenologically, the latter parameter con¬ 
figuration models the limiting case of a Sl-like scenario, 
where /i = 0. In such a dynamical process all individuals 
are infected in the steady state. 
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TABLE I. Structure and spectra of the normalized network 
of layers rj). The eigenvalues assumes that the average 

degree of each layer, (fc*), is the same, i.e. (fc*) = {k),Vl. 


Network 

4>J(A,77) 

Eigenvalues 

Line with 2 nodes 


1 + 

Line with 3 nodes 

Rfc'-i) 1 0 1 

L 0 2 

(k) 

(k) - 

(k) + 722 

Multiplex 

1 1 1 

{k)+n 


III. SPECTRAL ANALYSIS OF 7^(A, rj) 

As observed on the previous section, the supra adja¬ 
cency tensor TZ{X, rf) plays a major role on the epidemic 
process. Consequently, a deeper analysis of the spec¬ 
tral properties of such object can give us further insights 
about the whole process. First of all, the generalization of 
the eigenvector problem to the eigentensor is described 
on Appendix |A 2| allowing us to use some well estab¬ 
lished linear algebra tools. Additionally, in this section 
we generalize the spectral results of interlacing, obtained 
in [Ml 131, to the tensorial description adopted here. Be¬ 
sides, we also make use of the inverse participation ratio, 
IPR(A), as a measurement of eigenvalue localization |81| . 
As a convention, we assume that the eigenvalues are or¬ 
dered as Ai > A 2 > ...Anm and the individual layer eigen¬ 
values are denoted as A*. 


A. Interlacing properties 


Invoking the unique mapping presented on Ap¬ 
pendix A 2 and considering the results of [Ml 133] , we 
can use the interlacing properties to relate the spectra of 
the multilayer network with the spectra of the network 
of layers. First of all, we define the normalized network 
of layers in terms of the supra contact tensor as 




( 8 ) 


where we are implicitly assuming a multilayer network 
in which the layers have the same number of nodes and 
a dependency on the spreading rates (the demonstration 
that such tensor is an unfolding of the matrix exposed 
in |M] is shown on Appendix |A 3[ ). Additionally, let’s 
denote by > Mm the ordered eigenvalues of 

4)-(A, 77). Following |30| . the interlacing properties imply 


Anm. — 


nm—m-\ 


+3 — Mi S 


< A, 


(9) 


for j = TO,..., 1. As examples, Tablejljshows the spectrum 
of three simple networks of layers that can be computed 
analytically: a line with two and three nodes and a tri¬ 
angle. Figure 12 shows a schematic illustration of those 3 
multilayer networks. 

Furthermore, using similar arguments we can also ob¬ 
tain results for the normalized projection, formally given 
as 


(If) 

whose ordered eigenvalues, denoted by > 7/2 > ... > 
Vm., also interlace with the supra contact tensor satisfying 

Anm—n+j — <A„ (11) 

for j = n,..., 1. Finally, the adjacency tensor of an ex¬ 
tracted layer also interlaces, yielding 

Anm-n+j < A' < A^, (12) 

for j = n,..., 1. These results show that the eigenvalue of 
the multilayer adjacency tensor is always larger than or 
equal to all of the eigenvalues of the individual isolated 
layers as well as the network of layers. 

The interlacing properties presented here imply some 
constraints to the epidemic threshold. As advanced 
in [3Q], let Ai{Ai) be the i-th eigenvalue of the tensor 
Ai and consider that the set of eigenvalues is ordered as 
before. Moreover, for simplicity, we suppress the argu¬ 
ment when referring to the supra-contact matrix. First of 
all, assuming a fixed ratio of spreading rates, we observe 
that the eigenvalue of the multilayer follows 


/Ay^_ 1 ^ 1 

VmA Ai(A^(f)) - Ai’ 


Vfe 1,2 ,...,to, (13) 


where 



is the critical point for the single layer f and 


UJc " Ai( 4>|) - Ai’ 


(14) 


where denotes the critical point of the network of 

layers. Finally, considering the projection, we get 


UA "A^- A^’ 


(15) 


where is the critical point of the normalized pro¬ 

jection. Thus, the spreading process on the whole system 
is at least as efficient as it is on the layers and on the net¬ 
work of layers. Note that efficiency is understood here in 
terms of the position of the critical point, and not re¬ 
garding the fraction of infected individuals in the steady 
state. 
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FIG. 1. Schematic Illustration of the 3 multilayer networks cases considered as examples. Top panels represent the original 
networks which give rise to three distinct conhgurations for the networks of layers. See the text for more details. 


B. Localization and spreading of diseases 


Next, we investigate the behavior of the system near 
the phase transition and whether the phenomenon of dis¬ 
ease localization shows up. These two issues were ex¬ 
plored for monoplex networks in |in| and m , respec¬ 
tively, but have not been addressed for the case of multi¬ 
layer systems. The nodal probabilities can be written as 
a linear combination of the eigenbasis of TZ as 

A 

where c(A) are the projections of on the eigentensors 
/. Similarly to m, substituting such expression on the 
middle term of eq. |B7| we obtain 

^ ^ ^ AEA<c(A0A7ay(A0+A. ’ ^ ^ 


Considering only the contributions of the first eigen¬ 
value and eigentensor, for A > Ac, the first order approx¬ 
imation of the macro state parameter is p « cxit, where 
T = which yields 




(18) 


Such an expression is exact if there is a gap between the 
first two eigenvalues uniEi]. Furthermore, considering 
two eigentensors we have p « aiT + ■ Besides, fol¬ 

lowing a similar approach as in m we can use the inverse 


participation ratio: 

IPR(A)= (/^7A))V^^ (19) 

In the limit of nm oo, if the IPR(A) is of order 0{1) 
the eigentensor is localized and the components of fpg{A) 
are of order 0(1) only for a few nodes. On the other 
hand, if IPR(A) —)■ 0 then this state is delocalized and 

the components of /^|(A) ~ Additionally, 

another possible scenario, completely different from the 
traditional single layer one, is possible if we consider lo¬ 
calization on layers instead of on a fraction of nodes. In 
such a case, the IPR{A) will be of order 0{l/n) in the lo¬ 
calized phase, whereas it will be of order 0{llnm) in the 
delocalized phase. This is because, in the localized phase 
the components of the eigentensor are of order 0{l/y/ri) 
for all the nodes in the dominant layer and of order zero 
for nodes in other layers. Observing that, one easily real¬ 
izes that the correct finite-size scaling to take in order to 
characterize such a transition is m —>■ oo, i.e., the number 
of layers goes to infinity while the number of nodes per 
layer remains constant. In fact, in this limit IPR{A) will 
vanish on one side of the transition point while remain¬ 
ing finite on the other side. In this way, we can observe 
localized states also in the case in which there is no pos¬ 
sibility for localization in each of the layers if they were 
isolated. 


IV. MONTE CARLO SIMULATIONS 

We next compare the analytical results with Monte 
Carlo simulations of the spreading process. The method 
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proposed in mm is adapted here to the case of mul¬ 
tilayer networks. At each time step the time is incre¬ 
mented by At = ' (^jVi+AjVfc+> 7 jv^) '’ where N, is the num¬ 
ber of infected nodes, and Nk and Nm are the number of 
intra-layer and inter-layer edges emanating from them, 
respectively. With probability , one ran¬ 

domly chosen infected individual becomes susceptible. 
On the other hand, with probability 
infected individual, chosen with a probability propor¬ 
tional to its intra-layer degree, spreads the disease to an 
edge chosen uniformly random. Finally, with probabil¬ 
ity (^jv +Aivr+r;A'' ) infected individual, chosen with a 
probability proportional to its inter-layer degree, propa¬ 
gates the disease to an edge chosen uniformly. If an edge 
between two infected individuals is selected during the 
spreading, nothing happens, only time is incremented. 
The process is iterated following this set of rules, simu¬ 
lating the continuous process described by the SIS sce¬ 
nario. 

The quasi-stationary state (QS) method |111 155] re¬ 
stricts the dynamics to non-absorbing states. Every time 
the process tries to visit an absorbing state, it is substi¬ 
tuted by an active configuration previously visited and 
is stored on a list with M configurations, constantly up¬ 
dated. With a probability Pr a- random configuration on 
such a list is replaced by the actual configuration. In or¬ 
der to extract meaningful statistics from the quasi-static 
distribution, denoted by P{n^), where is the num¬ 
ber of infected individuals, the system must be on the 
stationary state and a large number of samples must be 
extracted. In this way we let the simulations run during 
a relaxation time tj. and extract the distribution P{n^) 
during a sampling time ta- The threshold can be esti¬ 
mated using the modified susceptibility na. given by 

where is the quasi-stationary distribution P{n^). As 
argued in [nisi] the susceptibility presents a peak at 
the phase transition on finite systems. Such measure is 
the coefficient of variation of the temporal distribution 
of states over time on the steady state. Note that the 
magnitude of the susceptibility x is not of primary inter¬ 
est to us, but rather the position of its maximum value 
with respect to /i/A, since it will coincide with the critical 
threshold for sufficiently large systems. 

In addition, after obtaining the curves of y x A by 
the QS method, we also apply a moving average filter in 
order to get rid of the noise. Such an approach improves 
the visual quality of the plots and does not interfere on 
the results, since the order of magnitude of the noise 
is smaller than those of the peaks corresponding to the 
transition points. 

The parameters used in the QS method are Pr = 0.01, 
ta varies from 10® to 10® and tr varies from 10® to 3 x 10® 
in order to obtain a smoother curve. The QS method 
demands a large sample size, since it is estimating the 


variance of a distribution. Moreover, we construct the 
X X A curves in steps of AA = 10“® and the moving 
average window has 5 points. 


V. 2-LAYER MULTIPLEX SYSTEMS 

In this section we numerically study 2-layer multiplex 
systems. First, we focus on the phase diagram of the 
spreading process as a function of the inter and intra 
layer spreading rates for both, SIS and SIR scenarios. 
Next, we analyze the spectral properties of such systems, 
comparing with results of Section |III[ Finally, we per¬ 
form Monte Carlo simulations that show the existence 
of multiple susceptibility peaks on multiplex networks. 
The latter results are analyzed in terms of the spectral 
properties of 7^(A, p). 


A. Numerical solution 


Results shown in this section are the numerical solu¬ 
tions of the ODE systems 1^ (SIS) and Cl (SIR) using a 
Runge-Kutta (4,5) algorithm |35|. We consider a 2 layer 
multiplex network (m = 2), where each layer has n = 10^ 
nodes. In order to build a multiplex network where the 
epidemic thresholds associated to the individual layers 
are well separated, we must guarantee that A)^ ^ A 2 - 
Therefore, we chose the degree distribution of the first 
layer to be P{k) ~ fc“^ ®, whereas that of the second 
layer is P{k) ~ Both layers are created using 

the uncorrelated configuration model m- Moreover, we 
consider a multilayer network in which every node has 
its counterpart on the other layer. This pairing of nodes 
of different layers is made randomly. Each result is the 
solution considering one single (and fixed) multiplex net¬ 
work. 

Figure[^shows the phase diagram considering the aver¬ 
age fraction of spreaders for the SIS dynamics (or recov¬ 
ered for the SIR dynamics) as the macro-state variable 
as a function of the spreading parameters A and rj for 
a given recovering rate ^ = 1. The dashed white line 
denotes the epidemic threshold obtained from eq. In 
(a) we show the SIS scenario, while (b) corresponds to 
the SIR model. In both cases, it is possible to observe 
two changes on the system’s behavior. The first on the 
epidemic threshold, while the second near the epidemic 
threshold of the second layer. In addition, we note the 
agreement between the theoretical epidemic thresholds 
and the numerical results. Furthermore, the higher rj, 
the lower the epidemic threshold, which is a consequence 
of the eigentensor problem. Also note that p increases 
for a fixed X as rj increases, even for A ~ 0, which means 
that in such extreme cases, the disease spreads mainly on 
the interlayer edges. 

Figure]^ shows the phase diagram for ^ = 1 and differ¬ 
ent values of the parameter rj for the SIS dynamics. For 
r] = 0 we have no inter-layer spreading, while for rj = 0.5 
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(a) SIS 





10 


(b) SIR 


1 0.48 
0.42 
] o .36 
0.30 
0.24 
0.18 
Jo.12 
0.06 
0.00 

1 0.72 
0.64 
0.56 
0.48 
0.40 
0.32 
0.24 
0.16 
0.08 
0.00 


FIG. 2. Phase diagrams over a 2-Layer multiplex system, 
where each layer is a scale-free network with n = 10 ^ nodes, 
for a fixed value of 71 = 1. (a) Density of spreaders as a 

function of the parameters rj and A. (b) Density of recovered 
individuals as a function of the parameters rj and A. Colors 
represent the fraction of spreaders and the white line is the 
threshold calculated using equation 


we have a fixed spreading rate, independent of the intra¬ 
layer rates. In addition, we also evaluated cases where 
the ratio ^ is constant. In Fig. |^(a) we have the global 
behavior of the system, which is an average of the indi¬ 
vidual behavior of the layers, represented in panels (b) 
and (c), since both layers have the same number of nodes. 
Furthermore, we also observe that the two individual net¬ 
works show different behaviors near the epidemic thresh¬ 
old [Tn|. The first layer (Fig.[^(b)) has a lower epidemic 
threshold than the second. However p grows (as a func¬ 
tion of A) slower than in the second. This feature can 
be observed clearly in Fig. (b) and (c), where we show 
results for 77 = 0 , that is, when there is no spreading 
between the layers. 

Considering the discrete system, Cozzo et al. ver¬ 
ified the shifting on the dominated layer (the largest 


amongst all individual eigenvalues) as the ratio ^ in¬ 
creases. Here we observe the same effect, as can be seen 
in Fig. (c). Additionally, we can also note another 
global change approximately beyond A > (A^)”^. Our 
findings suggest the possibility of multiple phase transi¬ 
tions due to the multiplex structure of the network. It 
is noteworthy that in spite of the similarities between 
our continuous model and the discrete model j^, both 
represent slightly different processes. On the continuous 
case, two events cannot happen at the same time. On the 
other hand, on the discrete model, every node contacts 
its neighbors on one discrete time step. Despite these dif¬ 
ferences, the results show that both the continuous and 
discrete formulations are phenomenologically similar. 


B. Spectral analysis 


Since the epidemic process is described through the 
supra adjacency tensor TZ{X,r]), its spectral properties 
give us some insights about the whole process, especially 
about the critical properties of the systems under anal¬ 
ysis. In this section we focus on the spectral analysis 
of such tensor as a function of the ratio j considering a 
2 -layer multiplex network with two different layers, i.e., 
there is a distance between the leading eigenvalues of each 
layer. Some important aspects of the spectral properties 
are left to Appendix where we present an analytical 
approach to the problem of eigenvalue crossings on Ap¬ 
pendix |DT)^ We focus on two special cases in increasing 
order of complexity: (i) the identical case, presented on 
Appendix |D 1 b| where both layers are exactly the same 
— i.e., there is a high correlation between the degree on 
each layer —; and (ii) the non-identical case, discussed in 
Appendix |D 1 c| where both layers have the same degree 
distribution, but different configurations. 

In this section we focus on the case of two different 
layer structures, with spaced leading eigenvalues. Con¬ 
sidering a multiplex network made up of two scale-free 
networks with 7 « 2.2 and 7 « 2.8. Both layers have 
{k) « 8 and n = 10 ^ nodes on each layer and the leading 
eigenvalues are A} = 42.64 for the first and A^ = 21.29 
for the second. 

Figure shows the spectral properties of the tensor 
TZ{X,r]) as a function of the ratio?. In contrast to the 
identical layers (see Appendix D 1 b) and the case of sta¬ 
tistically equivalent layers (Appendix Did, figures 10 
and [H where some eigenvalues increase while others de¬ 
crease, here all the observed eigenvalues always increase. 
Moreover, we do not observe any crossing or near-crossing 
behavior. Regarding IPR(A), the same pattern as for the 
similar case is found: for small values of j and consider¬ 
ing the first eigenvalue, the system appears localized on 
the first layer and delocalized on the second, while for 
IPR(A 2 ), it is the contrary. For larger values of both 
layers contribute equally to the IPR(A). Furthermore, 
the main difference we observe for the current setup with 


respect to the two similar networks (see Fig. 11 pre- 


















A A 

(a) Global (b) First layer 


8 



A 


(c) Second layer 


FIG. 3. Individual layer behavior over a 2-Layer multiplex system. Each layer has n = 10^ for a fixed value of /r = 1. The 
results considering both layers are shown in (a), while the dynamics in the individual layers are shown in (b) (P{k) ~ 
and (c) (P{k) ~ The arrows indicate the layers leading eigenvalues. 


sented on Appendix 0. is that now no drastic change 
on the inverse participation ratio is found, as expected, 
since there is no near-crossing. 

From figure we can also extract an important nu¬ 
merical result regarding the perturbation theory. We 
observed that in our case, considering a two spaced- 
individual layer eigenvalues problem, the leading eigen¬ 
value can be approximated by the largest leading eigen¬ 
value of the individual layers for y 1, such approx¬ 
imation becomes poorer as j increases, but it can be 
acceptable up to ^ 10, within a certain error. Apart 

from that, note that both eigenvalues tend to increase, 
while its difference tends to decrease. 

Furthermore, analyzing the eigenfunction properties. 
Fig. [5 shows the contribution of each layer to the IPR(A) 
considering different values of j. Results correspond 
to a multiplex network composed by two Erdos-Renyi 
networks, both with n = 5 x 10^, the first layer with 
(k) = 16, while the second has (k) = 12. Observe that 
for lower values of ^ the main contribution comes from 
one layer, configuring a localized state and consequently 
placed on one axis (the a;-axis) of Fig. Then, when 
the ratio ^ increases, there is a transition to a delocal¬ 
ized state. This corresponds to an increase of the inverse 
participation ratio of the second layer, however at the ex¬ 
pense of decreasing the value of the inverse participation 
ratio of the first layer. In other words, in the localized 
phase, only the entries of the eigenvector associated to 
the dominant layer are effectively populated, while the 
entries associated to other layers are not. In the delo¬ 
calized phase all the entries are equally populated. The 
inset of the figure further evidences this transition: it 
represents the angle, 9, between the vector composed by 
the IPR contributions, v = [lPR(Ai[), IPR(Af)]^, and 
the x-axis, where a change from zero to 45 degrees is ob¬ 
served as the ratio ^ is increased and the system goes 
from a localized to a delocalized state. 


C. Multiple susceptibility peaks 


Mata and Ferreira showed that it is possible to have 
multiple susceptibility peaks on monoplex networks |35| . 
They studied the behavior of a SIS model on networks 
with 7 > 3. Here we show that such phenomena also 
appear, in a natural way, on multilayer networks. Moti¬ 
vated by the findings reported in the latter sections, es¬ 
pecially by the presence of a second change in the slope 
of p as observed in figures and we have performed 
extensive Monte Carlo simulations using the QS-method 
with the aim of determining as accurately as possible the 
points at which the transitions takes place for a 2-layer 
multiplex network. Here we use the multiplex built up in 
Section V B[ since the leading eigenvalues of each layer 
are spaced. Note that our numerical simulations are per¬ 
formed on a fixed network, since we follow the quenched 
formalism. 


Figure]^ shows that for low values of the ratio j, both 
networks are weakly coupled and the system exhibits two 
well-defined susceptibility peaks (vertical dotted lines). 
However, as this ratio increases the peak signaling the 
presence of the second critical point decreases and even¬ 
tually vanishes. In our simulations, we have observed 
that up to j « 1, the second peak, although less defined, 
is still present. Beyond the latter point, only one peak re¬ 
mains. As ^ further increases, the position of the critical 
point remains the same, and the peak is even more well 
defined. Interestingly enough, if the ratio j continues 
to increase — in our case beyond y ^ 10 — the critical 
point shift to the left to values that are even smaller than 
the smallest critical point of the individual layers. It is 
worth highlighting that a similar qualitative behavior can 
be seen in the results shown in Fig. (a), where one can 
also observe a second change in the slope of p near the 
leading eigenvalue of the second layer. This change also 
vanishes as the intra-layer spreading increases. 

Since the tensor TZ{X,ri) plays a major role on the 
spreading process, our spectral results can help under- 
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FIG. 4. Spectral properties of the tensor TZ{X, rf) as a function 
of the ratio ^ for a multiplex with two layers, the first with 
7 « 2.2, while the second 7 « 2.8. Both have (fc) « 8 . On the 
top panel we present the inverse participation ratio (IPR(A)) 
of the two larger eigenvalues and the individual layer contribu¬ 
tions, while on bottom panel we show the leading eigenvalues. 
Every curve is composed by 10^ log spaced points, in order to 
have enough resolution. 


standing the observed critical dynamics. In epidemiolog¬ 
ical terms — or in general for contagion processes —, 
the localization of the disease on a certain layer means 
that most of the spreading is expected to take place on 
the nodes of that layer. Moreover, in addition to the lo¬ 
calization on the layers, one can also have localization 
effects on specific nodes or groups of nodes, for instance. 

In order to analytically explain this phenomenon, we 
evaluate IPR(A) for the two leading eigenvalues, as this 
measure indicates the localization of an eigenstate, see 
Section VB (results shown in Fig.|^. Comparing the sus¬ 
ceptibility and IPR(A), we observe that IPR(A 2 ) starts 
decaying for ? « 1 and crosses the value ^— , at which 



IPR{Al) 


10 ^ 

10^ 


10^ 

C'l'< 

10 ° 




10 '^ 

10 '^ 


FIG. 5. Diagram of the contribution of each layer to the 
IPR(A) for different values of the spreading ratio j. The 
dashed line represents the case where both layers have the 
same contribution, i.e. a line with slope one. In the in¬ 
set, we show the angle 6 between the vector composed 
by the contributions of each layer to the IPR(A), v = 
[IPR(A(), IPR(Ai)] and the x-axis. The multiplex network 
used here Is composed of two Erdos-Renyi networks, both with 
n = 5 X 10^, the first layer (k) = 16 ((A()~^ ~ 0.0625), while 
the second (fe) = 12 ((Ai)“^ « 0.0833). 


the associated eigenvector delocalizes, for ^ « 10, com¬ 
paring well with the point at which the second peak in the 
susceptibility decays and finally disappears. Moreover, 
IPR(Ai) decays from 3 < ^ < 10, which coincides with 
the range where the remaining maximum in the suscep¬ 
tibility reaches higher values and is better defined. More 
interestingly, note that IPR(Ai) is mainly composed by 
the contributions of the first layer for a lower spreading 
ratio, suggesting that it is localized on such layer. There¬ 
fore, our results suggest that the IPR(A) is a proper mea¬ 
sure to detect and predict the observed localization phe¬ 
nomena and potentially for m localization transitions, as 
we will show on Section IVIl 

Regarding the definition of a critical point it is im¬ 
portant to highlight that the concept of phase transition 
only applies in the infinite size limit (the thermodynamic 
limit). However, on the literature of complex network 
dynamics, specially for epidemic spreading, it is usual to 
use the terms critical point and phase transition on finite 
systems, since we find a behavioral change on that point. 
More importantly, for scale-free networks such point van¬ 
ishes in the thermodynamic limit. Following the usual 
convention on the complex network literature, the first 
susceptibility peak observed on all the experiments can 
be classified as a critical point of a phase transition. On 
such point, the dynamics goes from a disease-free state 
to an endemic state. On the other hand, the second sus¬ 
ceptibility peak cannot receive this classification, since 
the process is already on a endemic state. Although it 
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FIG. 6 . Susceptibility, x, as a function of the spreading rate A 
for different ratios of inter and intra-layer spreading ratings, 
j for a fixed value of /r = 1 over a 2-Layer multiplex system, 
where each layer have n = 10 ^, the first with 7 « 2 . 2 , while 
the second 7 « 2.8. Both have (k) ~ 8. The simulated values 
are 2 = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 
1.3, 1.4, 1.5, 1.6, 2, 3, 4, 5, 6 , 7, 8 , 9, 10, 20, 30. 


cannot be considered as a critical point, we have a tran¬ 
sition from a localized state to a delocalized state. In 
other words, before the second susceptibility peak most 
of the events take place on only one layer (the one with 
largest individual eigenvalue), while after this point both 
layers are active and spreading the disease. 


D. Second susceptibility peak analysis: 

Erdos-Renyi layers 

The second peak on the susceptibility curve suggests 
the existence of a second order phase transition. How¬ 
ever, from its existence alone we cannot conclude this 
unequivocally, since although this point is related to the 
delocalization of the disease, the system is already in an 
endemic phase (upper critical regime in Physics jargon). 
Observe that if r; € O (;^), in the thermodynamic limit 
we would have a phase transition. However, such con¬ 
figuration cannot be considered as a multilayer network, 
since both layers are (virtually) decoupled. Addition¬ 
ally, observe that we only analyzed layers without cor¬ 
relation. Such features can also introduce different phe- 
nomenologies, some were briefly explored in |32| . however 
for discrete-time. 

In order to better understand the second peak of sus¬ 
ceptibility we analyze a 2-Layer multiplex network com¬ 
posed by two Erdos-Renyi networks, in which we can 
precisely control the mean degree and consequently the 
epidemic threshold by fixing the number of edges. Fur¬ 
thermore, for scale-free networks with a divergent second 
moment of its degree distribution, the epidemic threshold 


vanishes in the thermodynamic limit [T]. On the other 
hand, Erdos-Renyi networks always have a non-zero and 
finite critical point. Aside from that, since the nodes 
on such a network are statistically equivalent, the proba¬ 
bilities are expected to be approximately the same. 
Henceforth we assume that the first layer has a higher 
connectivity, that is, a lower epidemic threshold. 

First of all, analyzing the layers individually for ^ > 
(A})”^ > A)"^ the first layer is in its upper critical regime 
(endemic state), while the second layer still is in its sub- 
critical regime (disease-free state). Then, for a coupling 
parameter, 77 > 0 , the probability of a node on the second 
layer being infected also increases. In fact, for Erdos- 
Renyi layers, it will be always larger than zero. There¬ 
fore, we can map this problem into an e-SIS model m, 
where each node has a probability of experiencing a spon¬ 
taneous infection. Note that such a model does not 
present an absorbing state. In this mapping, we are in¬ 
terested on the behavior of the second layer and consider 
that the self-infection e is determined by the contribution 
of the first layer by means of the contacts between nodes 
in different layers, which are Poisson processes with pa¬ 
rameter r]. This would imply that we would not have a 
second order phase transition. However, we have a transi¬ 
tion from a localized system, in which only the first layer 
is active and able to sustain the disease for long times, 
to a delocalized system, where both layers are active. 

In order to explore the time evolution of the system for 
a set of parameters near the second susceptibility peak, 
we run the continuous simulation 50 times and perform a 
moving average Alter over a sampling of the original time 
series, resulting in 5 x 10^ points. This approach give 
us an average curve over time. Note that for continuous 
simulations the number of points can vary from one run 
to another. Both networks used have n = 5x 10^, the first 
(k) = 16 ((Aj)“^ « 0.0625), while the second (k) = 12 
((Af)-i « 0.0833). 

Figure]^ shows the time evolution of a disease spread¬ 
ing on the second layer for different values of A and 77 . 
The initial conditions for these experiments consider that 
the first layer has an initial probability of a node being 
infected equal to 0 . 01 , while on the second every node 
is a spreader. Note that we chose this initial condition 
for visual purposes, since any initial condition would re¬ 
sult in a similar steady state regime. In this way, during 
the transient state we observe a decay of the fraction of 
infected individuals, then, at the meta state that config¬ 
ures the steady state, we observe a stochastic variation 
centered on the average value. Besides, such fluctuations 
tend to increase near a “critical point”. We observe that 
for (Aj)“^ ^ ^ ^ (^ 2 )”^ for 77 = 10 “"* the incidence 
is very low, of order O {^), however, larger than zero. 
As we increase the value of A we drive the system to its 
active state, being able to sustain the disease and spread¬ 
ing it by the intra-edges contacts. Besides, increasing 77 
we are able to increase the incidence of the disease due 
to the intra-edge contacts. Near the critical point of the 
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FIG. 7. Time evolution of the fraction of infected nodes on the second layer for p = 1, different values of ?; = 
10“^, 10~®, 10“^, 10“^ and different values of spreading rate: (a) A = 0.078, (b) A = 0.083, (c) A = 0.085 and (d) A = 0.088. 
The multiplex network used is composed of two Erdos-Renyi networks, both with n = 5 x 10'*, the first layer {k) = 16 
((A})“* ~ 0.0625), while the second (k) = 12 ((Af)“* ~ 0.0833). 


second layer, ^ = (A^) * = 0.833, we can observe some 
features that are similar to a transition. From below, we 
observe that the lower the value of rj, the longer it takes 
for the system to reach the steady state, similarly to what 
it is expected in phase transitions. On the other hand, 
slightly above the critical point, the time to get into the 
steady state decreases and the curves for rj = 10 “* and 
ry = 10“^ get closer. This suggests that the effects of 
intra-layer spreadings are the main source of spreading. 
Finally, for ^ sufficiently large, we observe the same be¬ 
havior for all values of 77, i.e. all of them are in an active 
state. 

In addition to the analysis shown in this section, we 
also inspected in detail the steady state for different sys¬ 
tem sizes, showing that neither the fluctuations diverge 
nor the final fraction of infected individuals goes to zero 
on the second layer. This analysis suggests that we do not 
have a second order phase transition but that the dynam¬ 
ics changes from a localized to a delocalized phase. In this 
reach phenomenological scenario, the transition point is 
still of great importance for practical purposes, for in¬ 
stance when it comes to study immunization policies. 
These complementary results are shown in Appendix|D 2| 


VI. 3-LAYER INTERCONNECTED SYSTEMS: 

THE BARRIER EFFECT 

Following the main ideas of the last sections, we ex¬ 
plore the spreading dynamics in multilayer networks with 
more than two layers. Specifically, we have carried out 
numerical simulations for a 3-layer system. We generate 
multiplex networks using three scale-free networks, with 
7 « 2.3, 7 « 2.6 and 7 « 2.9, with (k) « 8 and n = 10^ 
nodes on each layer. Note that we consider three lay¬ 
ers with spaced individual leading eigenvalues in order 
to investigate whether multiple susceptibility peaks are 


a generic phenomenon of multilayer systems. Note that 
we have two possible topologies for the network of lay¬ 
ers: (i) a line graph and (ii) a triangle (which is a node- 
aligned multiplex). In its turn, the first can be arranged 
in three possible configurations by changing the central 
layer. That is, we have four possible systems. In this sec¬ 
tion we focus on two configurations, the multiplex case 
and the line (2.3 -I- 2.9 -I- 2.6). Both cases summarize the 
richness of dynamical processes on interconnected net¬ 
works, presenting a new phenomenon, the barrier effect 
of an intermediate layer. We proceed by analyzing the 
spectral properties of this multilayer system in terms of 
the inverse participation ratio and the susceptibility. Re¬ 
garding the other interconnected networks, we present 
those complementary results and analyses in Appendix[E| 
Additionally, in Appendix |E 1| we show that increasing 
also increases the role of the inter-layer edges relative 
to the intra-layer ones. Consequently, the structure of 
the network of layers imposes itself more strongly on the 
eigenvalues of the entire interconnected structure. 


A. Spectral analysis 

Figure]^ shows the IPR(Ai) of tensor TZ. On the main 
panel we present the individual contribution of each layer, 
while on the insets we have the total IPR(Ai). On the 
top panel we have the line ( 2 .3-1-2.9-1-2. 6 ), whereas on the 
bottom panel we have the multiplex network. In this sec¬ 
tion we focus on the spectral comparison of two cases: (i) 
the lines (2.3-1-2.6-I-2.9) and (2.3-I-2.9-I-2.6) and (ii) the 
line (2.6-1-2.3-1-2. 6 ) and the multiplex network. Addition¬ 
ally, the reader is referred to Appendix |E 1| specifically 
to Fig. [^for complementary results. 

An interesting phenomenon can be observed comparing 
the different configurations of the network of layers. The 
largest eigenvalue of the whole system, Ai, has its associ- 
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FIG. 8. Spectral properties of the tensor TZ{X, r]) as a func¬ 
tion of the ratio ^ for a multiplex with two layers with the 
same degree distribution (different random realizations of the 
configuration model) and connected to its counterpart on the 
other layer. On the top panel we present the inverse partic¬ 
ipation ratio (IPR(A)) of the two larger eigenvalues and the 
individual layer contributions, while on bottom panel we show 
the leading eigenvalues. Every curve is composed by 10® log 
spaced points, in order to have enough resolution. On (a) we 
have the line (2.3 -I- 2.9 -I- 2.6), while on (b) the multiplex case. 


ated eigenvector localized in the dominant layer, that is, 
in the layer generated nsing 7 = 2.3. Regarding the line 
configuration, depending on the position of that layer 
in the whole system — i.e., central or peripheral layer 
— the contribution of the non-dominant layers to the 
IPR(Ai) varies. In particular, when the dominant layer 
corresponds to an extreme node of the network of layers, 
the contribution of the other two layers will be ordered 
according to the distance to the dominant one. Conse- 


FIG. 9. Susceptibility y as a function of A considering all 
three layer configurations and many different ratios which 
is represented by the color of the lines. The recovering rate 
is ^ = 1. The simulated values are ^ = 0.05, 0.06, 0.07, 0.08, 
0.09, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 2, 3, 4, 5, 6, 
7, 8, 9, 10, 20. On (a) we have the line (2.3-I-2.9-|-2.6), while 
on (b) the multiplex case. 


quently, when the dominant layer is in the center of the 
network of layers, the contributions of the non-dominant 
ones are comparable (see Fig. 17 on Appendix E1 for 
complementary results). 

Furthermore, for the first eigenvalue, which is usually 
enough to analyze the localization as a first order ap¬ 
proximation, we observe that the layer with the largest 
eigenvalue dominates the dynamics. In addition, note the 
similarities between the multiplex and the line configu¬ 
ration (2.6 -I- 2.3 -b 2.6) (see also Fig. 17 Appendix Ell, 
where the non-dominant layers behave similarly. This 
is because for small values of the effect of the extra 
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edge in the network of layers (closing the triangle) is of 
order 77 ^ and so the similar behavior observed for the 
two configurations. As ^ grows, the symmetry in the 
node-aligned multiplex dominates the eigenvector struc¬ 
ture and the contributions of all layers are comparable. 
As we next show, the different contributions of the lay¬ 
ers to the total IPR(Ai) are at the root of the multiple 
susceptibility peaks observed. 


B. Multiple susceptibility peaks 


Figure shows the susceptibility as a function of A for 
different ratios of y. We observe three well defined peaks 
on such curves when the ratio j is small. In addition, 
similar to the 2 -layer case, such peaks tend to become 
less defined and vanish as the ratio j increases. The 
third peak is less defined than the others because the 
average number of infected nodes is larger in this case. 
Consequently the susceptibility tend to be lower, since it 
measures the variance in relation to the average. Such an 
observation suggests that it could be harder to observe 
peaks for non-dominating layers that have an individual 
critical point too far from the dominating layer. 

Except for the line (2.3-1-2.9-1-2. 6 ) all figures are similar 
and present similar peaks, implying that the susceptibil¬ 
ity peaks occur approximately at the same point (for a 
complementary analysis see Appendix E3 and Fig. 18 1 . 
On the other hand, the line (2.3 -I- 2.9 -I- 2.6) shows a 
slightly different behavior for the second peak, that is 
found for a larger value of A than for the other cases. 
This result suggests that when the layer with the largest 
eigenvalue is located at the center of the line, it can ef¬ 
fectively act as a barrier to the disease. In addition, it is 
verified that the extra inter-edges of the multiplex case 
does not lead to radical changes on the transition points. 
We remark that the susceptibility does not measure the 
fraction of spreaders in the steady state. Thus, despite 
of the similarities of those curves, the phase diagrams for 
the incidence of the disease are different. 

Coming back to what is observed for the network of 
layers described by the line (2.3 -I- 2.9 -I- 2.6), an inter¬ 
esting phenomenon arises, namely, the formation of bar¬ 
riers to the epidemic spreading. Since the middle layer 
has the lowest individual eigenvalue among the layers, it 
creates a barrier effect “delaying” the second transition. 
Moreover, we observe that this transition also vanishes 
for higher values of the ratio y, if compared to the other 
cases. This can be related to the inverse participation 
ratio of Ai, IPR(Ai), shown in Fig.|^ Note that, for the 
line ( 2 .3-1-2.9-1-2. 6 ), the contribution of the layer 7 = 2.6 
is the lowest. As shown in Section [V A| (and in |28]), for a 
2 -layer multiplex, the non-dominant layer has its critical 
point shifted to a lower value of the spreading rate, which 
means that the outbreak takes place before it would have 
happened if that layer were isolated. However, here such 
shifting is compromised by the fact that the central layer 
is unable to sustain the epidemic process, acting effec¬ 


tively as a barrier for disease contagion. Apart from this 
new effect, the system behaves qualitatively similar to 
the 2 -layer scenario. 


VII. CONCLUSIONS 

In this paper, we have generalized and extended previ¬ 
ous analyses to the case of multilayer networks. To this 
end, we have made use of the tensorial representation 
introduced in |29| . which allows to extract upper and 
lower bounds for the disease incidence of a SIS model 
and the critical points for both, the SIS and the SIR dy¬ 
namical processes. We have also validated our analytical 
insights with extensive numerical simulations, recovering 
results like those presented in [55] regarding the shift¬ 
ing of the global epidemic threshold to lower values of 
the spreading rate and the role of the so-called dominant 
layer. Furthermore, we have observed a transition on 
the spectra of the supra-contact tensor, from the spectra 
resulting from the union of the individual layers to the 
spectra of the network of layers. This behavior implies 
that other dynamics and more complex structures can 
also be significantly affected by the interconnected na¬ 
ture of the system. In addition, we have also character¬ 
ized analytically the phenomenon of eigenvalue crossing 
on the supra-contact tensor for the case of two identical 
layers. It is worth noticing that any dynamical process 
that is described by the same matrix will be affected by 
this effect. 

Our main results concern the emergence and vanishing 
of multiple susceptibility peaks as a function of the ratio 
between the inter-layer and intra-layer spreading rates 
and their relation to the spectral properties of the mul¬ 
tilayer, which also revealed the phenomenon of disease 
localization, and in particular, its relation with the exis¬ 
tence of crossings or near-crossings of eigenvalues. Using 
the QS-Method and Monte Carlo simulations, we have 
been able to precisely determine the transition points. 
We remark that the first susceptibility peak is a phase 
transition, from a disease free state to an endemic, but 
localized, state. On the other hand, the second peak is a 
transition from a localized to a delocalized state, which 
is not a second order phase transition. Additionally, we 
have proposed an analytical approach based on the use of 
the inverse participation ratio to characterize such transi¬ 
tions as a localization phenomenon, thus also connecting 
with |31| . 

A detailed exploration of the parameter space showed 
that as the ratio between the inter-layer and intra-layer 
spreading rates increases, the peaks of the susceptibility 
measured for the non-dominant layers tend to occur at 
lower values of A and vanish as ^ increases up to a point 
in which only one susceptibility peak is observed, which 
is a true phase transition. Interesting enough, our results 
point out that such a transition can take place for even 
lower values of A than the inverse of the largest leading 
eigenvalue among all individual layers. 
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Finally, another important finding presented here is 
the opposite phenomenon, namely, the barrier effect, 
which happens when the susceptibility peak takes place 
at a larger value of A than that expected as a consequence 
of the multiplex topology. Specifically, if the layers are 
arranged in such a way that the one with the smallest 
leading eigenvalue is at the center of the network of layers 
(for instance, as it happens for the line (2.3 + 2.9 + 2.6) 
configuration), then the corresponding transition could 
be delayed due to the barrier effect. Summarizing, our 
results emphasize the importance of studying multilayer 
systems as they are and not only as a collection of indi¬ 
vidual layers. 
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Appendix A: Tensorial representation 

In this appendix we extend some important concepts of 
the tensorial representation. On Section ] A l| we present 
the projections, while on Section |A 2 j we show the equiv¬ 
alence of the eigentensorial problem and the eigenvec¬ 
tor problem of the supra adjacency matrix. Finally on 
Section |A 3| we prove the relation between the tensorial 
projection and the matricial representation, which is fun¬ 
damental to the interlacing results. 


where T/ S Note that such a network presents 

0 

self-loops, which are weighted by the number of edges on 
the layer. Additionally, since we assume that the layers 
have the same number of nodes, the edges of the network 
of layers have weights equal to the number of nodes n. 

Another important reduction of the multilayer network 
is the so-called projection [53]. Such network aggregates 
all the information into one layer, including self-loops 
that stand for the number of layers in which a node ap¬ 
pears. Mathematically, we have 

= M^^Ul (A2) 

where e 

2. Eigenvalue problem 

As presented on the main text, the epidemic threshold 
is closely related to the leading eigenvalues of the supra- 
contact tensor. Here we describe the eigenvalue problem 
considering the tensorial representation. Such eigenvalue 
problem can be generalized to the case of a rank-4 tensor 
leading to 

7^;|/ay(A)=A/^^(A), (A3) 

where A is an eigenvalue and is the correspond¬ 

ing eigentensor. In addition, we are assuming that the 
eigentensors form an orthonormal basis. Importantly, 
the supra-contact matrix, R, in |28| can be understood 
as a flattened version of the tensor TZ“J{X,r]). Gonse- 
quently, all the results for R also apply to the tensor 
TZ. As argued in |29| . that supra-adjacency matrix cor¬ 
responds to unique unfolding of the fourth-order tensor 
m yielding square matrices. Following this unique map¬ 
ping we have the correspondence of the eigensystems. 
Here, we consider that the eigenvalues are ordered as 
Ai > A 2 > ...Anm and denote the individual layer eigen¬ 
values as A(. 

3. Proof of Equation |8| 


1. Tensorial projections 

For the sake of completeness we present other projec¬ 
tions of multilayer networks, which are specially conve¬ 
nient on tensorial notation, due to its compactness. Be¬ 
sides the adjacency tensor presented on the main text, 
the network of layers |3D] also characterizes the topology 
of the system. In this reduced network representation, 
each node represents one layer and the edges between 
them codify the number of edges connecting those two 
layers. Formally we have, 

'kl = , (Al) 


Gonsidering the matricial representation of a multi¬ 
layer network, given by 


A = ©aA“ + C = 


Al C12 

C21 A2 


Clm 

C2m 


Cml C, 


m2 


An 


(A4) 


where A G A“ G is the adjacency matrix 

of the layer a G { 1 , 2 , ...m} and C is a coupling matrix. 
Since we assume multilayer network in which the lay¬ 
ers have the same number of nodes we have Cij = I. 
Assuming a partition of such network, represented by 
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S G which is the characteristic matrix of such 

partition, where Sij = 1 if i G Vj and zero otherwise, 
where Vj is the network of layers partition. 

In order to use the results of 123 [M] we have to prove 
that the network of layers matrix H isniisii is an unfold¬ 
ing of our tensor ‘I’lCA, 77 ), formally given by 

E = r-^S'^AS, (A5) 


where F is a diagonal matrix with normalizing constants 
(for more, see references [SOI IM]b In words, the prod¬ 
uct AS is a summation over the blocks of the matrix A, 
resulting in a matrix with the degree of each node. The 
subsequent left product with impose another summa¬ 
tion, whose result is a matrix composed by the sum of 
all elements of the blocks. Finally, the product by F”^ 
normalize the result by A Formally we have. 




. 

• fcl™ ■ 





^5 = 

A,ml 

km2 . 

jxmm 


(A6) 


where G is a vector with the number of edges 

emanating from each node on layer i to layer j and AS G 
Then, 






EA:22 . 


x;a:™i 




(A7) 


where ^ G K are scalars with the number of edges 
that connect a node on layer A to a node on layer j. Fi¬ 
nally, the product by F“^ introduce the average degree 
instead of the summation, producing the same results as 

Eq.i 


Appendix B: The Susceptible-Infected-Susceptible 
(SIS) model analysis 

In this section we present an extension of the analysis 
presented on the main text regarding the SIS model. To 
begin with, we present comments on the exact model def¬ 
inition and its relation with the first order approximation 
on Section [BT| On Section [B^ we present a derivation of 
the critical point for the first order approximation, while 
on Section |B 3| we present the derivation of the lower and 
upper bound for such model. 


1. Model definition: complementary comments 

In probability theory and stochastic processes it is 
usual to define random variables as capital letter. How¬ 
ever, that is the same usual notation for tensors. In order 


to avoid confusion we will use bold capital letters for ran¬ 
dom variables. For instance, we define the Bernoulli ran¬ 
dom variable that defines the state of a node as S^g, 
where it assumes one of two values, zero if the node 
(36 is susceptible or one if it is infected. By definition, 
Xps = iSps), where (•) is the expectation operator and 
X^g is the probability of the node (35 being infected. 

In this way, without any assumption on the indepen¬ 
dence of random variables the exact equation can be writ¬ 
ten as 


+ (1 - Sfgs) , 

(Bl) 

where the supra contact tensor is defined in|^ This equa¬ 
tion can be interpreted as an exact version of the epi¬ 
demic process m- However without any approximation 
the solution of such problem involves 0 ( 2 "™) equations, 
since we have to write the expressions for the expecta¬ 
tion for all the products. The first order approximation 
consists in (S^gSa^) « {Spg){Sa^) = Such 

approximation is shown on eq. Interestingly, observe 
that eq. |Bl| can be written in terms of the covariance, de¬ 
fined as Cov[S'^^, S'ay] = {SpgSa^) - {S^g){Sa^). Con¬ 
sequently, isolating the probability of the product and 
substituting it in eq. |B 1 | we find, by inspection, that the 
error is given by Cov[S'^^, Say], which is assumed to be 
zero. In |39| , the authors observed this relation and pro¬ 
posed an accuracy criteria for monoplex networks. 


2. The epidemic threshold 


An important concept for dynamical systems that 
present an absorbing state and an active phase is the 
critical point. Considering the SIS process, below this 
point the system is inactive and the disease tends to dis¬ 
appear. On the other hand, for above this point we have 
the active phase, where the disease is present on a frac¬ 
tion of the population. Assuming /i > 0 and that the 


dynamics has reached the steady state, 
can write eq. as 


dt 


= 0 , we 


l-X 





^oc'-y voo 


(B2) 


Expanding the left-hand term following the geometrical 
series, where < 1 , we 

obtain 

00 ^ 

(DE(^S) (B3) 

k^l 

In addition, similarly to [H], suppose = ef^g, 
where e is an arbitrary small constant and f^g > 0. Sub¬ 
stituting in eq. |B3| and dividing by e we have 

7^;|(A,r7)/ay = (0/^^ + e(0 ' + O(e^). (B4) 
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Considering a sufficiently small e > 0 this expression re¬ 
duces to the eigentensor equation 

= (B5) 

leading to the critical point 

( 9 ^ = a.. m 

where Ai is the largest eigenvalue of TZ, which is the same 
as the largest eigenvalue of R in |28j . 


which can be inserted into Eq. |B7| to give, 
> 1 - ^ 


1 -f 


(a) 


A \ ^min _ I 


(B12) 


Finally, combining Eqs. |B8| and |B12| the bounds of 
Eq. are 


1 - 


1 + ^ 


(J) 


- 1 




(^) dps + 1 


(B13) 


3. Upper and lower bounds for the steady-state 


In order to obtain some bounds for the epidemic inci- 

dX, 


dence considering the steady state, where 


'-/35 


dt 


— 0. For 


a monolayer system those bounds were calculated in m- 
We consider a multilayer network without self loops and 
denote the steady state of each node as X'^. Then, im- 


dX 


posing 


VQO _ 


ps 


dt 

A7^“|(A,77)X° 

A7^;|(A,7y)X“- 


S' 


= 0 to Eq. 3 we have 


= 1 - 


A7^;J(A,r;)X-, + l 


(B7) 


The value of X'^- is then obtained by iterating the above 
equation from an initial value, until convergence. Upper 
and lower bounds can be obtained by considering only 
the first iteration of Eq. |B7| For the upper bound we 
have 




1 


{t) + 1 


(B8) 


where 


dp~s = 'i^;]i\v)Uc.p = 




(B9) 


As can be noticed, there are two different contribu¬ 
tions to the upper bound coming from intra and inter¬ 
layers connectivity. Both of them tend to increase the 
probability of a node being infected. Furthermore, the 
higher is the degree, the higher is this upper bound. 
On the other hand, for the lower bound, let us denote 


Min{X-j} = X” 

we have 


Then, substituting X™™ in Eq. 


B7 


AT” 


> 1 - 




(A,?7)U„.^X™-kl 


(BIO) 


Denoting Minjc?^^} = d™”, we obtain 

1 


AT" 


> 1 - 


(^) 


(Bll) 


Appendix C: The Susceptible-Infected-Recovered 
(SIR) Model 

Aside from the SIS epidemic model, we can also con¬ 
sider the SIR model. Contrasting with the SIS, which 
have just one absorbing state (inactive), the SIR have 
many absorbing states. In fact, considering an infi¬ 
nite population we have an infinite number of absorbing 
states. 


1. Model definition 


Introducing the recovered and susceptible states, here 
denoted by Ygj and ^/3(5i respectively. Then, using a sim¬ 
ilar notation as in the latter section and associating Pois¬ 
son processes to nodes and edges, we have the dynamical 
set of equations 


d^ps 

dt 

d^ps 

dt 

d^ps 

dt 


-dXps + Zps^n;]{x,v)Xo.p 

d’Z^pS 

-ZpsXn‘;]{x,v)Xc.p. 


(Cl) 


Note that the Poisson processes on the nodes model the 
recovering, whereas on the edges, model the spreading. 


2. Epidemic threshold 


Since there is no dynamic steady state in the SIR 
model, the epidemic threshold has a different interpre¬ 
tation from that of the SIS model. Above the threshold 
the total number of recovered individuals reaches a fi¬ 
nite fraction of the population, when the dynamic starts 
with a small fraction of infected individuals. Formally, 
the initial condition are: Ar^^(O) = = 0 and 

Zpsi^) ~ ^ ~ 7771 where c is a small constant, c <C nm. 
Neglecting higher order terms, we have 


dAT, 


ps 


-MA^^,~ + A7^;|(A,r?)A„^. 


(C2) 


dt 
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After a proper factorization, 



where (5“J is a tensor analogous to the identity matrix, 
whose elements are one if the indices are the same. The 
epidemic threshold is as in eq. which is the critical 
value for both SIR and SIS dynamics. 


Appendix D: 2-Layer Multiplex systems 

In this section we present some complementary analy¬ 
sis for the 2-Layer multiplex case. Here we focus on some 
spectral aspects of such systems, mainly on the eigen¬ 
value crossing and near-crossing phenomenon, presented 
on Section |D 1| and additionally on the second suscepti¬ 
bility peak using a finite size analysis. Such results are 
presented on Section |D 2| and are complementary to Sec¬ 
tion on the main text. 


1. Spectral aspects 


In this section we focus on the spectral analysis of 
the tensor TZ{X,r]) as a function of the ratio First 
of all, we present an analytical approach to the problem 
of eigenvalue crossings on Section |D 1 a| then we focus 
on three special cases in increasing order of complexity: 
(i) the identical case, where both layers are exactly the 
same. Thus, there is a high correlation between the de¬ 
gree on each layer, presented on Section D 1 b[ (ii) the 
non-identical case, where both layers present the same 
degree distribution, but different configurations on Sec¬ 
tion |D 1 c| The case of two different layer structures, 
considering that their leading eigenvalues are spaced was 
presented on the main text. 


a. Eigenvalue crossing 

Let us analyze the spectra of a simple setup: multiplex 
networks composed by I identical layers. Such class of 
networks provides insights about the spectral behavior as 
a function of (^). Although they are not very realistic a 
priori, there are situations in which this representation is 
helpful: for instance, in the context of disease contagion, 
one might think of a multi-strain disease in which each 
strain propagates in a different layer allowing co-infection 
of the host population. 

The adjacency tensor can be written as 




V 

A 


FIG. 10. Spectral properties of the tensor TZ{X, rf) as a func¬ 
tion of the ratio ^ for a multiplex with two layers with the ex¬ 
act same degree distribution and connected to its counterpart 
on the other layer. On the top panel we present the inverse 
participation ratio (IPR(A)) of the three larger eigenvalues, 
while on bottom panel we show the leading eigenvalues. Ev¬ 
ery curve is composed by 10^ log spaced points, in order to 
have enough resolution. 


Kronecker delta. Observe that the sum of two Kronecker 
products, A = Ira® A + ® In, where is the iden¬ 

tity matrix of size n and Km is the adjacency matrix of 
the complete graph with m nodes is the unfolding of the 
adjacency tensor in this case. In this way, the eigenvalue 
problem can be written as 


7^;|(A,r;)=A^<5| + |J^A|, (Dl) 

where is the 2-rank layer adjacency tensor, is 
the adjacency tensor of the network of layers, which is 
a complete graph on the multiplex case, and 5^ is the 


(D2) 

where the sum of the eigenvalues of A, A\, and K, /i^, are 
also eigenvalues of the adjacency tensor, hence = 
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b. Identical layers 


Considering a multiplex network made up of two layers 
with the same configuration. Each layer of the multiplex 
is a network composed by n = 1000, (k) « 6, A* = 14.34, 
with degree distribution P{k) ~ k~^ ’^. Aside from the 
intra-edge configuration, we also impose that inter-edges 
connect a node with its counterpart on the other layer, 
i.e., every node has the same intra-degree on all layers. 
Such a constraint imposes a high correlation between the 
degrees on each layer. 

Figure shows the spectral behavior of such a mul¬ 
tiplex as a function of the parameter (2^. On the top 
panel, we represent the inverse participation ratio of the 
first three eigenvalues, while on the bottom panel, we 
plot the first ten eigenvalues. When the ratio ^ = 0 the 
eigenvalues have multiplicity two, as can be seen on the 
left side of the bottom panel (approximately, since the 
figure starts from 10“^). More importantly, those eigen¬ 
values tend to behave differently: one increases, while 
the other tends to decrease. This behavior leads to the 


eigenvalue crossing (see Appendix D 1 a). The inset of 


the bottom panel zooms out the region where the cross¬ 
ing takes place. Note that the eigenvalues cross at the 
same value for which the inverse participation ratio shows 
an abrupt change. Indeed, the jump in the IPR(A) has 
its roots in the interchange of the eigenvectors associated 
to each of the eigenvalues that are crossing. Moreover, 
we stress that the abrupt change observed for IPR(A) is 
always present in such scenarios, but it could be either 
from the lower to the higher values or vice versa depend¬ 
ing on the structure of the layers. 


FIG. 11. Spectral properties of the tensor IZ{X,r]) as a func¬ 
tion of the ratio ^ for a multiplex with two layers with the 
same degree distribution (different random realizations of the 
configuration model) and connected to its counterpart on the 
other layer. On the top panel we present the inverse partic¬ 
ipation ratio (IPR(A)) of the two larger eigenvalues and the 
individual layer contributions, while on bottom panel we show 
the leading eigenvalues. Every curve is composed by 10® log 
spaced points, in order to have enough resolution. 


(A( + y/ij) faj, i = 1,2, ...n and j = 1, 2, ...to. Then, 

(A( + ^M,) = (Ai + ^/A.). (D3) 

The eigenvalues of the complete graph are /ii = to — 1, 
and yti = —1, Vz > 1, yielding to 


V K-K 

X TO 


(D4) 


which imposes crossings on the eigenvalues of the adja¬ 
cency tensor for identical layers, since (^) is a continuous 
parameter. 


c. Similar layers 

In addition to the identical case, we have also consid¬ 
ered a multiplex network composed by two layers with the 
same degree distribution (i.e. the same degree sequence), 
with P{k) ^ but different random realizations of 

the configuration model. Furthermore, the inter-edges 
follow the same rule as before, connecting nodes with 
their counterparts on the other layer assuring that ev¬ 
ery node has the same intra-degree on all layers. Each 
layer of the multiplex network is composed by n = 1000 
and (k) « 6. Since each layer is a different realization of 
the configuration model, both present a slightly different 
leading eigenvalue, the first Aj = 15.21 and the second 
Af = 14.34. 

Figure shows the spectral behavior of such a mul¬ 
tiplex in terms of the largest eigenvalues, on the bottom 
panel, and the IPR(A), on the top panel. Here, in ad¬ 
dition to the global inverse participation ratio, we also 
present the contribution of each layer to this measure. 
Such analysis is meaningless on the identical case, since 
the contribution is the same. As shown in the figure, 
we observe that for small values of in regard to the 
first eigenvalue, the system is localized on the first layer 
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FIG. 12. Final number of infected nodes on the second layer (with lowest individual eigenvalue) as a function of the size of the 
layers on the main panels, while on the insets we present the fraction of infected nodes on the left and the standard deviation 
on the steady state on the right. The parameters used on the simulations are shown on the tile of each panel. They are a 
combination of the parameters A = 0.078,0.083,0.085,0.088 and rj = lO”"^, 10“®, 10“^, 10“^. Furthermore, the layer sizes are 
n = 2 X 10^3 X 10^4 x 10^ 5 x 10^ 6 x 10^ 7 x 10^ 8 x 10^9 x 10^10^,2 x 10'*, 3 x 10'‘,4 x 10* and 5 x 10* and m = 2 on 
all cases. 


and delocalized on the second. On the other hand, the 
picture changes when we focus on the second eigenvalue, 
as it is localized on the second layer, but delocalized on 
the first. For larger values of j, both layers contribute 
equally to IPR(A). Analogously to the identical case. 


there is a change on IPR(A 2 ), which seems to be related 
to the changes on A 2 , as one can see on the bottom panel 
and in the inset. Note that for this case, there is no cross¬ 
ing, i.e., the eigenvalues avoid the crossing -also referred 

































































































20 



n 


n 


FIG. 13. Final fraction of infected nodes on the layer with 
lowest individual eigenvalue as a function of the the size of 
the layers. The colors represent different values of 77 , while 
on we have A = 0.078 on (a), A = 0.083 on (b) A = 0.085 
on (c) and A = 0.088 on (d). Furthermore, the layer sizes are 
n = 2 X 10^ 3 X 10^ 4 X 10^ 5 X 10^ 6 X 10^ 7 x 10^ 8 x 10^ 9 x 
10^ 10'‘, 2 X 10^, 3 X lO'*, 4 X 10'‘ and 5 x lO'* and m = 2 on all 
cases. Each curve is the result of a parameter rj, from bottom 
to top ri = 10 "^, 10 ■^ 10 "^ 10 “^ 


to as near-crossing. 


2. Finite size analysis 

In this section we analyze the behavior of a 2-layer 
multiplex network at the steady state considering differ¬ 
ent sizes. Such a multiplex was built considering two 
Erdos - Renyi networks with a fixed mean degree. As 
mentioned in the main text, we chose this type of net¬ 
works because their epidemic threshold do not vanish at 
the thermodynamic limit, which contrasts with the scale- 
free networks. In this way, we have a well-defined critical 
point that can be precisely tuned regardless of the net¬ 
work size. Following the usual convention on the complex 
network literature, the first susceptibility peak observed 
on our experiments can be classified as a critical point 
of a phase transition. On such point, the dynamics goes 
from a disease-free state to an endemic state. However, 
the second susceptibility peak cannot be classified as a 
second order phase transition, since the disease is already 
in an endemic state. Although it cannot be considered 
as a critical point, before the second susceptibility peak 
most of the events take place on only one layer (the one 
with the largest individual eigenvalue), while after this 
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FIG. 14. Finite size analysis of the susceptibility. On the 
main panel we have the susceptibility as a function of A for 
different sizes of 2 layer multiplex network, where the first 
layer have (k) = 16 and the second (k) = 12. On this ex¬ 
periment we fixed the ratio ^ = 0.01. On the inset we show 
the susceptibility of the two peaks as a function of the layer 
size, where the blue symbols refer to the Hrst peaks, while the 
green symbols refer to the second peak. Besides, the red lines 
are a linear fitting of those points. The layer sizes evaluated 
are n = 3 X 10^ 4 X 10^ 5 x 10^ 6 x 10^ 7 x 10^ 8 x 10^ 9 x 
10^ lO'*, 2 X 10‘‘, 3 X 10"*, 4 X 10'*, 5 x 10*, 10®. 


point both layers are active and spreading the disease. 


Similarly to the experiments shown in Section VD 


here we run the continuous simulation 50 times and per¬ 
form a moving average filter over a sampling of the orig¬ 
inal time series, resulting in 5 x 10* points. The sim¬ 
ulations are run up to t = 10^. Note that for contin¬ 
uous simulations the number of points can vary from 
one run to another. The steady state statistics are es¬ 
timated for t > 950 or in other words, the last 50 time 
units. In contrast with the main text, here we are in¬ 
terested in comparing results for different network sizes, 
n = 2x 10^, 3 X 10^, 4 x 10^, 5 x 10^, 6 x 10^, 7 x 10^, 8 x 
10^, 9 X 10^, 10*, 2 X 10*, 3 X 10*, 4 x 10* and 5 x 10* and 
m = 2 in all cases. Besides, we considered the mean de¬ 
gree as (k) = 16 for the first layer and (k) = 12 for the 
second. We expect that the second susceptibility peak 
appears near the epidemic threshold of the second layer 
individually, i.e. A « 0.083. 


Figure [^presents the number of infected nodes in the 
steady-state on the layer with the lowest individual eigen¬ 
value as a function of the size of the layers and a com¬ 
bination of the parameters A = 0.078,0.083,0.085,0.088 
(near the individual critical point of the second layer) 
and 77 = 10“*, 10“^, 10“^, 10“*. Besides, on the insets 
we have the information about the average fraction (left 
inset on each panel) and its fluctuations, measured by 
the standard deviation (right inset on each panel). The 
straight lines in red were obtained by a least squares re- 
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FIG. 15. Distribution of the eigenvalues. On the rows, from top to bottom, for the interconnected networks of Lines 2.3+2.6+2.9, 
2.3 + 2.9 + 2.6, 2.6 + 2.3 + 2.9 and the multiplex. On the columns, from left to right, we varied the ratios ^ = 1,10,100 and 
1000 respectively. All histograms were built with 100 bins. 


gression method. 

We observe an approximately linear behavior of the 
number of infected nodes on the second layer as a func¬ 
tion of the number of nodes on such layer (see the main 
panels of Fig. 12 1 . Consequently, the fraction p 2 also 


presents a linear trend (see the left inset on each panel 
of Fig. 12). In fact, it presents a flat pattern, i.e approx¬ 


imately constant. Besides, the number of infected nodes 
is always larger than zero, since it is not a disease-free 
state. Furthermore, we also observed that the fluctua¬ 
tions tend to be very low (see the right inset on each 
panel of Fig. 12). Regarding the fluctuations, it is note¬ 


worthy that on a phase transition they tend to diverge, 
which does not happen in our analysis, thus also rul¬ 
ing out a second order phase transition as far as it con¬ 
cerns. We also note that fluctuations are slightly higher 
for lower spreading rates, as can be seen by the error bars 
for A = 0.078, which is explained by the delocalization of 


our system. 

Furthermore, in figure we present the comparison 
of steady state fractions. In each panel, we fix a value of 
A and compare different values of t]. It emphasizes the 
influence of rj on the final fraction of infected nodes on 
the second layer. Note that for rj = 10“^ and small net¬ 
works the behavior exhibits a growing trend. This is due 
to the fact that for networks with n < 10 ^ the contribu¬ 
tion of the first layer can be effectively neglected. In fact, 
observe that for n > 10 ^ the fraction of infected nodes 
on the second layer follows a flat pattern (see Fig. (c) 
and (d)). Finally, in figure 14 we present a finite size 


analysis of the susceptibility for different sizes, ranging 
from n = 3 X 10^ to n = 10® and m = 2 layers. Each 
curve was obtained using the QS algorithm, with which 
we simulated 120 points from A = 10“^ to A = 10“^. 
On such experiment we fixed the ratio j = 0.01. Addi¬ 
tionally, we also used a moving average Alter with two 
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FIG. 16. Evaluation of the 8 first eigenvalues of TZ{X, rf) for 
the multiplex configuration as a function of of the ratio 
It is noteworthy that such plot is visually equivalent for all 
the layer topologies composed by 3 layers. The dashed lines 
represents the individual layer leading eigenvalues. 


points for visualization purposes. In the inset, we show 
the scaling of the susceptibility corresponding to the two 
peaks. The positive slope for the first peak indicates that 
it divergences as the system size goes to infinity, thus ev¬ 
idencing the phase transition. On the other hand, the 
curve for the case of the second peak is flat whatever the 
value of the system size is, indicating that in contrast to 
the behavior observed for the first peak, in this case there 
is no divergence in the thermodynamic limit nor the peak 
vanishes. 


Appendix E: 3-Layer interconnected systems: 
complementary analysis 


In this section we study the introduction of a third 
layer, which increases the complexity of the system al¬ 
lowing four different network layer configurations, the 
line, which has three different configurations depending 
on the position of the layers, and the triangle, which is 
also a multiplex. This section is organized as follows: in 
the first subsection we perform the spectral analysis of 
the adjacency tensor as a function of the parameter 
showing that as we increase this parameter the spectral 
distribution tends to the spectra of the network of lay¬ 
ers, which is explained by interlacing theorems. Next, on 
sections |E 2| and |E 3| we show the complementary results 
of localization and susceptibility analysis, respectively. 
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FIG. 17. Spectral properties of the tensor 77.(A, 77 ) as a func¬ 
tion of the ratio ^ for a multiplex with two layers with the 
same degree distribution (different random realizations of the 
configuration model) and connected to its counterpart on the 
other layer. On the top panel we present the inverse partic¬ 
ipation ratio (IPR(A)) of the two larger eigenvalues and the 
individual layer contributions, while on bottom panel we show 
the leading eigenvalues. Every curve is composed by 10^ log 
spaced points, in order to have enough resolution. 


1. Spectral analysis 


Since the epidemic process is described through the 
supra adjacency tensor TZ{X,ri), its spectral properties 
give us some insights about the whole process, especially 
about the critical properties of the systems under analy¬ 
sis. Moreover, as the structure of the network of layers is 
not trivial anymore, we shall find important differences 
regarding the spectra of such tensors for the different 
topologies of the network of layers. 
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FIG. 18. Susceptibility x as a function of A considering all 
three layer configurations and many different ratios j, which 
is represented by the color of the lines. The recovering rate 
is ^ = 1. The simulated values are j — 0.05, 0.06, 0.07, 0.08, 
0.09, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 2, 3, 4, 5, 6, 
7, 8, 9, 10, 20. 


Figure shows the spectrum of the four configura¬ 
tions of networks when varying the ratio ^ = 1,10,100 
and 1000. Observe that we do not show the ratio y = 0 
since it is just the union of the individual layers’ spec¬ 
trum. For ^ = 1, the four configurations are very similar, 
especially the line graphs. In such case, the inter-layer 
edges are treated in the same way as the intra-layer ones. 
In other words, they are ignored and the network can be 
interpreted as a monoplex network. As the spreading 
ratio increases the spectrum tends to be clustered near 
the values of the eigenvalues of the network of layers. 
Such spectra was analytically calculated in Section [II| and 
shown in Table U on the main text. 


Regarding the triangle configuration, the clustering of 
the spectrum as j increases is even clear. Triangles 
present the lowest eigenvalue with multiplicity two. On 
the extreme case of ^ ^ see Fig. 15 we have 2/3 of the 


values near the left extreme value vmile 1/3 is near the 
leading eigenvalue. On the other hand, for the line config¬ 
urations, the frequencies of the eigenvalues distribution 
is related to the position of the central layer. However, 
on the limiting cases such differences are reduced. This 
pattern is naturally related to the increase of the spread¬ 
ing ratio: When ^ increases, so does the role of the inter¬ 
layer edges relative to the intra-layer ones. Consequently, 
the structure of the network of layers imposes itself more 
strongly on the eigenvalues of the entire interconnected 
structure. This comes as a consequence of the interlacing 
theorems shown in Section [III Al on the main text. 

Our findings can be related to the structural tran¬ 
sition shown in m, where the authors evaluated the 
supra-Laplacian matrix as a function of the inter-layer 
weights. Their main result is an abrupt structural tran¬ 
sition from a decoupled regime, where the layers seem 
to be independent, to a coupled regime where the lay¬ 
ers behave as one single system. Here, we are interested 
in the supra-adjacency tensor, however, we found a simi¬ 
lar phenomenological behavior and a structural change of 
the system as a function of the inter-layer weights, which 
in our case are determined by a dynamical process. 


2. Localization on interconnected networks 


Complementary to the results presented in Section [VT| 
here we present results for the lines (2.3 -I- 2.6 -I- 2.6) and 
(2.6 -I- 2.3 -I- 2.6). Similarly, the experiments here are 
conducted in terms of the inverse participation ratio, as 
it was done for the 2-Layer multiplex case. 

Figure shows the 10th larger eigenvalues of the 3- 
layer multiplex case. The dashed lines represent the lead¬ 
ing eigenvalue of each layer. Note that the leading eigen¬ 
value of the layer with P{k) ~ ® is the 7th larger 

on the network spectrum when ^ = 0. We observe that 
there is no crossings on the observed eigenvalues, which is 
an expected result, since the layers have different struc¬ 
tures. Furthermore, it is important to remark that all 
networks of layers evaluated also show similar qualitative 
behaviors. The topology of the network of layers does 
not lead to qualitative differences on the dependence of 
Ai on ^ for the first ten eigenvalues. We also notice that 
although it is only an approximation, the perturbation 
theory would be valid roughly up to ^ ^ 10. 

Figures 1^ and 17 shows the IPR(Ai). On the main 
panel we present the individual contribution of each layer, 
while on the insets we have the total IPR(Ai). As men¬ 
tioned on the main text, the first eigenvalue is usually 
enough to analyze the localization as a first order approx¬ 
imation. Here we observe that the layer with the largest 
eigenvalue dominates the dynamics. In addition, note the 
similarities between the multiplex and the line configu- 
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ration (2.6 + 2.3 + 2.6), where the non-dominant layers 
behave similarly. This is because for small values of 
the effect of the extra edge in the network of layers (clos¬ 
ing the triangle) is of order 77 ^ and so the similar behavior 
observed comparing the panel (b) of figures]^ and 17 for 


the two configurations. As ^ grows, the symmetry in the 
node-aligned multiplex dominates the eigenvector struc¬ 
ture and the contributions of all layers are comparable. 
As we next show, the different contributions of the lay¬ 
ers to the total IPR(Ai) are at the root of the multiple 
susceptibility peaks observed. 

Complementing and reinforcing the analysis of Sec¬ 
tion |VT1 comparing the different line configurations of 
the network of layers, observe that the largest eigenvalue 
of the whole system, Ai, has its associated eigenvector 
localized in the dominant layer, that is, in the layer gen¬ 
erated using 7 = 2.3. Depending on the position of that 
layer in the whole system — i.e., central or peripheral 
layer —, the contribution of the non-dominant layers to 
IPR(Ai) varies. In particular, when the dominant layer 
corresponds to an extreme node of the network of lay¬ 
ers, the contribution of the other two layers will ordered 
according to the distance to the dominant one. Conse¬ 
quently, when the dominant layer is in the center of the 
network of layers, the contributions of the non-dominant 


ones are comparable -note that in panel (b) of Fig. 17 


there is no difference in the contribution to IPR(Ai) of 
layers generated using 7 = 2.6 and 7 = 2.9. 


3. Multiple susceptibility peaks: additional results 

Figure [ 1 ^ shows the susceptibility as a function of A for 
different ratios of As observed in the main text, we 
also have three well-defined peaks in these curves when 
the ratio j is small. In addition, similar to the 2-layer 
case, such peaks tend to become less defined and vanish 
as the ratio ^ increases. 

Regarding the third peak, note that it is less defined 
than the others because the average number of infected 
nodes is larger in this case. Consequently the suscepti¬ 
bility tends to be lower, since it measures the variance in 
relation to the average. The comparison of Figures [^(b) 
and 18 shows that there is no difference in the position of 
the susceptibility peaks. As mentioned in the main text, 
the only observed difference is the barrier effect, shown 
in Fig. (a). We also remark the similarities between 
the line (2.6 -I- 2.3 -I- 2.6) and the multiplex case, which 
emphasize the role of the central node. In that line con¬ 
figuration, the layer with 7 = 2.3 spreads its influence 
to both layers, being this similar to the multiplex case, 
however with less intra-edges. 
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